function check_core_period(CONSTS, plot_data)

    c     = CONSTS.c;
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    d     = CONSTS.d;

    zeta_step = 0.05; % 0.05
    zeta_max = 1;
    zeta_vec = (zeta_step : zeta_step : zeta_max)'*d; 
    integ_m = zeros(size(zeta_vec));
    m_vec = (101:2:103)';
    for i = 1:size(zeta_vec, 1)
        zeta = zeta_vec(i);

        for im = 1:size(m_vec, 1)
            m = m_vec(im);

            period = pi/(k0*a);
            dots_per_period = 20; % 100
            order = 3;
            q_high = 0.5*log(10)*order/(k0*zeta);
            q_vec = (0.0 : period/dots_per_period : q_high)';

            y_vec_1 = func_1(q_vec, zeta, m, CONSTS);
            y_vec_2 = func_2(q_vec, zeta, m, CONSTS);
            integ_1 = trapz(q_vec, y_vec_1);
            integ_2 = trapz(q_vec, y_vec_2);

            integ_m(i, im) = - ((2*pi*k0^2*a)/c) * integ_1 ...
                           + (2*pi/(c*a)) * (1/m^2) * integ_2;

        end
    end
    
    if (plot_data)
        figure; hold on;
        integral_legend = {};
        m_sz = size(m_vec, 1);
        for im = 1:m_sz
            component = 0.5*(m_sz-im+1)/m_sz;
            plot(zeta_vec/d, integ_m(:, im), '.-', ...
                'Color', [component component min([component+0.5 1])]);
            integral_legend = [integral_legend; sprintf('Numerical integral m=%d', m_vec(im))];
        end
        hold off;

        title('K_m(\zeta)');
        xlabel('\zeta /d'); ylabel('K_m'); ylim([min(min(integ_m))*1.1 0]);
        legend(integral_legend);
    end

end